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SECTION  1  -  INTRODUCTION,  HEAT-FLUX  DATA  REDUCTION 


This  report  documents  the  data  reduction  code  referred  to  as  ACQ.  The  code  was 
originally  developed  at  the  Massachusetts  Institute  of  Technology  (MIT)  in  1984  for  the 
United  States  Air  Force  to  reduce  data  from  MIT’s  (then  recently  developed)  double-sided 
heat-flux  gauges.  The  ACQ  code  has  been  used  at  MIT  for  neaily  ten  years. 

The  study  conducted  at  Calspan  to  document  ACQ  has  several  different  facets.  The 
first  facet  examines  the  theory  and  methodology  applied  in  the  data  reduction  code.  The 
second  facet  of  the  study  ascertains  the  code's  cqiabUities  and  accuracy.  These  first  two 
facets  of  the  study  yielded  an  understanding  of  the  code  which  made  it  possible  to  apply  the 
code  in  a  more  user  friendly  and  versatile  way.  The  last  facet  of  the  study  documents  this 
application  of  the  code  within  the  Data  Acquisition  System  (DAS). 

The  ACQ  reduction  code  is  capable  of  reducing  data  from  both  single-sided 
(semi-infinite)  and  double-sided  (finite)  heat-flux  gauges.  Single-sided  gauges,  sometimes 
referred  to  as  thin-film  semi-infinite  gauges,  were  in  use  prior  to  the  development  of  the 
double-sided  gauge  at  MIT  and  are  still  being  used  at  many  laboratories  with  success. 

The  semi-infinite  gauges  consist  of  a  thin-film  resistance  thermometer  bonded  to  an 
insulating  substrate.  It  is  assumed  that  the  substrate  appears  infinitely  thick  to  thermal 
waves  propagating  from  the  surface.  Therefore,  the  time-dependent  heat-flux  can  be 
inferred  using  a  one  dimensional  heat  conduction  model.  The  semi-infinite  gauges  are 
typically  constructed  of  platinum  painted  onto  a  sufficiently  thick  piece  of  Pyrex  glass 
which  is  then  inserted  into  the  metal  airfoil.  Although,  in  some  instances  the  entire  airfoil  is 
constructed  from  the  insulating  material;  usually  a  machinable  ceramic. 
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The  double-sided  gauges  are  a  relative  of  the  semi-infinite  gauges  for  which  a 
second  thin-film  resistance  thermometer  is  bonded  on  the  opposite  side  of  the  insulating 
substrate.  This  eliminates  the  need  to  assume  that  the  substrate  is  semi-infinite  because  the 
measured  bottom  temperature  is  the  second  boundary  condition.  Therefore,  the  substrate 
ran  be  made  very  thin  making  it  possible  to  cover  the  airfoil  surface  with  a  layer  of  the 
substrate  material,  A  schematic  cross  section  of  a  double-sided  heat-flux  gauge  is  given  in 
Figure  1.1. 
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Figure  1.1  --  Schematic  Cross  Section  of  Double-Sided  Gauge 


The  double-sided  gauges  act  as  a  thermal  shunt  below  a  certain  fi-equency  where  the 
temperature  difference  across  the  insulator  is  a  direct  measure  of  the  heat-flux.  The  shunt 
cut-off  frequency  increases  as  the  insulator  thickness  is  reduced.  Conversely,  above 
another  frequency  the  substrate  appears  semi-infinite  to  the  upper  sensor  and  a  one 
dimensional  assumption  can  be  used  to  infer  the  heat-flux.  These  limiting  cases  are  not  in 
themselves  applicable  in  the  reduction  of  data  because  most  of  the  signal  frequencies  lie 
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between  the  limiting  frequencies.  Therefore,  a  numerical  technique  is  used  to  reconstruct 

the  entire  frequency  domain  front  DC  to  100  kHz. 

In  the  numerical  technique,  the  time  unsteady  one  dimensional  heat  conduction 
through  the  gauge  substrate  is  modeled  as  current  flow  through  a  discrete  resistance- 
capacitance  network.  A  schematic  of  this  analogue  is  given  in  Figure  1.2.  The  current 
through  the  first  resistor  represents  the  surface  heat-flux  for  the  measured  surface 
temperatures  which  are  the  input  voltages  to  the  circuit.  The  resistor  and  capacitor  elements 
are  analogous  to  the  thermal  properties  of  the  substrate.  The  circuit  elements  are 
logarithmically  distributed  through  the  network  to  match  the  insulator  attenuation 
characteristics  at  high  frequencies  to  increase  the  numerical  efficiency. 
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The  theory  applied  in  the  data  reduction  code  is  detailed  in  Section  2.  The  electric 
analogue  is  developed  with  comparison  to  an  alternative  finite-difference  approach.  The  R- 
C  circuit  constmction  is  outlined  along  with  the  solution  procedure  for  the  resulting  system 
of  equations. 

Section  3  documents  the  verification  of  the  ACQ  code.  MIT's  original  verification 
study  is  duplicated  along  with  additional  studies  to  ascertain  the  accuracy  of  the  code.  The 
two  test  cases  used  by  MIT  demonstrated  the  ability  of  ACQ  to  produce  heat-flux  data  from 
double-sided  gauge  data.  A  third  test  case  is  given  to  demonstrate  how  ACQ  can  be  used  to 
reduce  single-sided  semi-infinite  gauge  data.  The  second  part  of  the  verification  process 
examines  how  anomalies  in  the  recorded  data  affect  the  code’s  solution.  Truncation  error 
due  to  low  A/D  resolution  and  random  signal  noise  are  both  investigated. 

Section  4  outlines  the  application  of  the  data  reduction  code.  ACQ  has  been 
modified  to  enable  input  and  output  of  Data  Acquisition  System  (DAS)  format  files.  The 
program  can  be  operated  in  two  modes:  interactive  and  automatic.  The  interactive  mode 
provides  the  user  with  a  convenient  tool  for  generating  heat-flux  time  histories  for  one  or 
two  gauges.  The  automatic  mode  is  more  like  a  batch  process  which  not  only  performs 
multiple  gauge  heat-flux  calculations,  but  also  places  pertinent  information  into  the  run  log 
files.  Proper  application  of  the  code  is  addressed  which  requires  some  consideration  in  the 
design  of  the  experiment. 
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SECTION  2  -  DATA  REDUCTION  THEORY 


The  heat-flux  time  history  is  deduced  from  the  recorded  gauge  surface 
temperatures.  The  gauge  surface  temperatures  are  used  as  boundary  conditions  in  solving 
the  diffusion  equation  through  the  substrate  between  the  temperature  sensors.  The  thermal 
diffusion  equation  is  given  by  Equation  2.1  and  the  corresponding  heat-flux  is  given  by 
Equation  2.2. 


8T^  k  8^T 
^  pCp  dx^ 


(2.1) 

(2.2) 


2.1  —  Electric  Analogue 


An  electric  analogue  to  the  thermal  diffusion  equation  is  used  by  Norton  [1985]  in 
the  data  reduction  code  ACQ  supplied  to  Calspan  by  MIT.  This  analogue  is  convenient 
because  it  follows  the  theory  developed  by  Burd  and  Doe  [1980]  under  the  direction  of 
Oldfield.  The  predecessor  of  the  electric  analogue  circuit  described  by  Burd  and  Doe  was 
developed  by  Skinner  [I960].  In  this  theory  the  gauge  substrate  is  analogous  to  a 
resistance-capacitance  circuit  with  the  diffusion  equation  taking  the  form 


dV  ^  1  d^V 
W  rc  0x2 


(2.3) 


where, r=resistance/unit length, c=capacitance/unit length,  V  =  T,  Isq,cs pCp,  r  =  ^ 
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2.1.1  --  Discrete  Analogue 


A  discrete  R-C  circuit,  shown  in  Figure  2.1,  is  constructed  to  discretize  the  spatial 
domain.  For  the  semi-infinite  case,  Rn+1  is  set  to  a  very  large  value. 
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Figure  2.1  —  R-C  Circuit  Discrete  Analogue 


Consider  the  element  of  this  circuit  shown  in  Figure  2.2. 
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Figure  2.2  —  R-C  Circuit  Element 


From  elementary  circuit  analysis. 


Vi_i-Vi  =  IiRi 


Vi-Vi,i  =  l2Ri,i 

■3-,f 
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I3  =  li-l2 


dVj  Vi_,-V, 

*  dt  -  Ri  Ri,i 


rearranging. 


<IV|  V|-|  V|M  ,  I  ' 
<ft  Rfi  ^Ri  Ri,iJ 


V. 


i±i 


(2.4) 


Equation  2.4  is  a  first  order  Ordinary  Differential  Equation  (ODE)  for  the  voltage  at 
each  node  along  the  circuit.  To  demonstrate  that  Equation  2.4  is  truly  analogous,  an 


identical  ODE  is  derived  from  a  finite  difference  analysis  of  the  thermal  diffusion  equation. 
Consider  the  temperature  to  be  known  at  discrete  spatial  locations  which  are  not  equally 
spaced.  The  second  order  spatial  derivative  of  the  temperature  is  approximated  by 


a^Tj^  Ti_t  f  1  ,  IK 

3x2  AXilAXi  Axi+J  AXj+iAXj 


where, 


Axi  =  Xi-Xi_i 

^i+l  =  Xi+i-Xi 


substituting  into  Equation  2.1, 
Ti-i 


3Ti_ 


Ti 


(pC^O^  pCpAx.^ 


AX; 


1  + 


4+1 


i+1 


Ax 


i+1 


(pCpAxi) 


(2.5) 


Equation  2.5  is  analogous  to  Equation  2.4  with  the  capacitance  equivalent  to  the 
total  thermal  capacity  of  a  depth  Axj,  and  the  resistance  equivalent  to  the  total  thermal 
impedance  of  a  depth  Axj . 


Cj  =  pCpAXj 
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The  discrete  R-C  circuit  analysis  results  in  exactly  the  same  equation  as  the 
finite  difference  analysis  of  the  thermal  diffusion  equation.  Thus,  we  are  solving  the 
correct  analogous  equation.  Equation  2.4. 


2.1.2  -  Matrix  Theory  for  System  of  Equations 


The  application  of  Equation  2.4  to  all  the  nodes  of  the  circuit  produces  a  set  of 
coupled  first  order  ODE’s.  The  voltages  at  the  ends  of  the  circuit  are  input  and  analogous 
to  the  substrate  surface  temperatures.  Applying  Equation  2.4  at  the  endpoints, 

dVi  _  Vq  ViM  1  \  V2 
dt  RiCi  C^Rj  R2J  R2C1 

_  ^n-1  _  Ynl-L.  +  1  ^ 

dt  ^nl^n  ^n+li  ^n+l^n 


Casting  the  set  of  equations  into  vector  form. 


—  =  [A]V  + 
dt  RjC, 


Vo  r 
— +- 


_n±l 


^n+l^n 


where. 


(2.6) 


Vo  is  analogous  to  the  upper  surface  temperature  and  Vn+1  is  analogous  to  the  lower 
surface  temperature. 
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J_f  JL+_L 

Cj  Rj 


C.R, 

ifl 


^2V^2  ^3 


C2R3 


0  •  • 


[A]  = 


0  0  • 
0  0  0  0 


.J_(  J_+ J_ 

C.  1  R„  R. . 


This  coefficient  matrix  is  cast  into  a  symmetric  matrix  using  the  transformation 

y.  -  -Iji...  ( 

■ 


Equation  2.4  becomes 


dUj^  Uj-i  1  ^ 

dt  Ri/CAT  CilRi  Ri+ii 


^i+lV^^iQ+1 


The  set  of  equations  in  matrix  form  become 


^=[B]u+  ''""r-K 

dt  R,Vq 


where  [B]  is  a  tri-diagonal  matrix  with  the  following  elements 


^U+l  ^i+U' 


Rj+iv/CA+i 


Note  that  Vq  and  Vn+i  are  not  transformed  to  U  values  because  Cq  and  Cn-t-l  are 


undefined. 
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The  set  of  ODE's  are  uncoupled  by  defining  a  new  set  of  coordinates,  y,  such 


U  =  [Q]y 

where  [Q]  is  the  matrix  of  eigenvectors  for  [B].  Substituting  into  Equation  2.8, 


Multiply  through  by  [Q]‘ 


[Qr‘[Q]f =lQr'WQ]y+[Qr^£.  +[Qr^^^k.  am) 


V.O  i 

^n+lV^ 


0  •  •  0  ■ 

0  ^2  0  •  0 

Now  [Q]  *[®][Q]=  0  0  0  0 

0  0  0 

_0  •  •  0  V 

where  the  diagonal  elements  are  the  eigenvalues  of  [B].  Since  [B]  is  symmetric,  its 
eigenvectors  are  orthogonal,  hence 


Qll  Q21  Q31 

•  •  qm 

qi2  * 

qn2 

qn3 

[Qr=[Qr  = 


L^ln  Q2n  ^3n  *  *  ’  ^nn  J 

Then,  if  yi  is  the  i^h  element  of  column  vector  y.  Equation  2.10  gives, 


+  ^=r-+Qn.in 

Kiv't-i  i^n+lV'-'n 


(2.11) 


Thus,  Equation  2.1 1  is  the  uncoupled  ODE  to  be  solved  at  each  of  the  circuit  nodes. 
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2.2  Solution  Procedure 


The  theory  outlined  up  to  this  point  has  required  specific  resistance  and  capacitance 
values.  Since  the  thermal  properties  and  thicknesses  will  vary  from  gauge  to  gauge,  so  will 
the  specific  resistance  and  capacitance  values.  Thus  to  avoid  evaluating  the  eignevalues  and 
eigenvectors  for  each  gauge,  a  normalized  R-C  circuit  is  set  up  with 

(2.12) 


Then  Equation  2.11  becomes 


dyj^  1 
dt 


Vn.l  ^  1 
n+lV^i 


(2.13) 


and  the  transformation  given  by  Equation  2.7  becomes 


(2.14) 


To  have  the  same  eigenvalues  and  eigenvectors,  each  gauge  analogue  in  the  series 
must  have  the  same  number  of  circuit  nodes  (stages)  and  the  same  geometric  variation  in 
the  circuit  components.  The  specific  thermal  and  physical  properties  of  each  gauge  are 
introduced  through  Ri  and  Cl. 


2.2.1  —  R-C  Component  Selection 

Consider  the  circuit  shown  in  Figure  2.1  to  consist  of  a  series  of  'T'  sections  as 
shown  in  Figure  2.3. 


11 


Ri+2,1 

■sAAA — r 


I  Ri,1 

r-WNA 


Ri.2 

W\A 


Ri+1,1 

■WNA 


Ri+1,2 

W\A 


Ci 


Ci+1 


Ci+2 


Figure  2.3  -  R-C  Circuit  as  a  Series  of  T  Sections 


Let  P  be  the  resistor  matching  coefficient 


Ri,l  +  f^i,2 

(2.15) 

The  ratio  R/C  for  each  T  section  must  be  equal  to  r/c,  therefore 

r  ^i,l  +  ^i.2 

Ci 

(2.16a) 

+  (2,,6b) 

Ci+1 

The  geometric  variation  of  the  capacitance  and  resistance  is  chosen  to  be  logarithmic 

to  simulate  the  attenuation  characteristics  of  the  gauge  substrate. 

Ci+1  ~  ^Ci 

(2.17) 

Ri,2  f^i+1.1  ^i+1.2 

Ri.r  Ru  Ri.1.1 

(2.18) 

From  Equations  2.16  and  2.17, 

Rj+l.l  ^i+1.2 

J^i.l  +  I^i.2 

Then,  with  Equation  2.18  and  2.19 

0 

- =  Y  .•.  y=a^ 

1  +a 


(2.19) 


(2.20) 
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Replacing  Ri,2  in  Equation  2.16a,  we  obtain  the  relationship  for  the  first 
resistance  in  the 'T'  section. 

Ru=(§)tt!^  (2.21) 

Combining  Equations  2.15, 2.16a  and  2.21  gives 

(2-22) 

The  relationship  for  the  second  resistance  in  the  T'  section  can  be  written  as 


Returning  to  the  original 
be  defined  as  follows. 


(2.23) 


'L'  section  circuit  shown  in  Figure  2.1,  the  resistances  can 


(2.24) 

(2.25) 

(2.26) 


2.2.1 .1  -  Tx)garithmic  Geometric  Factor 

The  logarithmic  geometric  factor,  y,  is  determined  through  a  required  bandwidth 
consideration.  The  Bode  plot  for  a  finite  analogue  is  given  in  Figure  2.4. 
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Into 

Figure  2.4  -  Finite  Analogue  Bode  Plot 


In  the  normal  operating  range  an  analogue's  impedance  should  ideally  satisfy 


At  high  ftequencies  the  first  capacitor  shorts,  so 

|z»|  =  R, 

At  low  frequencies  the  impedance  of  the  capacitors  swamps  the  resistors,  so 


where  Q  is  the  total  capacitance  of  the  analogue.  The  useful  bandwidth  is  taken  as  the 
intersections  of  these  three  curves. 


Therefore, 


(2.27) 


\ 
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Defining  the  forcing  function  and  evaluating  the  bandwidth,  we  get 


pp_  ^max 


®min 


then 


FF  =  -^ 
c  Ri 


From  the  resistance-capacitance  relationship  defined  earlier  in  Equation  2.24 

FF  =  (1  +  Vy)^  (2.28) 

The  total  capacitance  is  the  sum  of  the  individual  capacitances  which  are  given  in  relation  to 
Cl  by  Equation  2.17. 

C,  =  ^C,  =  ^C,  +  IC,  +  fC,  +  ■•  ■  +  y"-‘C, 

i=l 

C,=C,^1  +  Y+y'+”+7”-' 


c,  r-i 

"C,  Y-1 


Substituting  into  Equation  2.28  yields, 

^^(l  +  Vy)(Y  -l)  ^2.29) 

Y-1 

Equation  2.29  is  solved  iteratively  for  the  geometric  factor,  y ,  with  the  desired  bandwidth 
and  specified  number  of  stages,  n. 

The  maximum  frequency  is  taken  as  one  half  the  sample  frequency  for  both  the 
semi-infinite  and  the  finite  cases.  The  minimum  frequency  is  taken  to  be  the  frequency  at 
which  the  thermal  wave  propagation  depth  equals  the  sensor  thickness.  This  minimum 
frequency  is  determined  according  to  the  following  analysis. 


e 


—  for  X  =  d 
e 
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where  d  is  the  gauge  substrate  thickness.  Solving  for  the  frequency, 


®min  = 


rkV  1 


(2.30) 


u;  pCpk 

Since  the  object  is  to  set  up  a  normalized  analogue,  the  pCpk  and  k/d  values  used  are  those 
for  the  first  sensor  to  be  analyzed. 


2.2.2  -  Setting  up  the  Normalized  Circuit 

The  normalized  resistance  and  capacitance  values  used  are 

R  j  =  1  and  C  j  =  1 

c;=7c;_i 

R;=VY(l  +  VT)c'i-i 

The  normalized  total  resistance  and  capacitance  are  calculated  as 


R 


m 

'.=Er. 


i=l 


c;=X'=‘ 


i=l 


(2.31) 

(2.32) 

(2.33) 


(2.34) 


(2.35) 


where  m=n  for  semi-infinite  and  n-i-1  for  finite.  These  totals  are  required  to  determine  Ri 
and  Cl  for  each  individual  gauge. 

The  normalized  [B]  matrix  is  constructed  and  stored  as  two  column  matrices  [D] 
and  [E].  The  diagonal  elements  of  [B]  are  contained  in  [D], 


f 


1  1 
4-- 


and  the  off-diagonal  elements  are  contained  in  [E], 
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The  eigenvectors  and  eigenvalues  of  [B]  are  calculated  by  a  subroutine  which  has 
[D],  [E]  and  an  identity  matrix  [I]  passed  to  it.  Upon  return,  the  [D]  column  matrix 
contains  the  eignevalues  and  the  identity  matrix  contains  the  eigenvectors. 


2.2.3  -  Dimensional  Solution 

The  physical  properties  of  each  gauge  are  introduced  to  determine  the  heat-flux 
corresponding  to  the  surface  temperature  history  of  each  particular  gauge. 


2.2 .3.1  —  Gauge  Physical  Properties 

The  physical  properties  are  introduced  through  Ri  and  Ci  for  each  gauge.  For  a 
finite  analogue,  the  analysis  yields  the  following  relationships: 

AX: 


by  definition 


R,= 

■  k 

q^pCpAx, 


d  1 

R.  = - r 

‘  k  R, 


^  d  pCpk 
‘  k  C 


(2.36) 


(2.37) 


For  the  serai-infinite  analogue  R)  and  C)  are  set  to  give  the  desired  maximum 
frequency.  Recall  Equation  2.27, 
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where 


Recall  Equation  2.24, 


(0 


max 


r  1 

cRl 


r _ 1_ 

c  pCpk 


'  VpCpka)o>ax 


R  —  — — 

^■c(i+V7) 


(i+Vy)V^ 


(2.38) 


(2.39) 


2.2.3.2  -  Integration  of  ODE 

The  solution  is  started  by  setting  the  initial  condition  at  each  node  in  the  circuit.  The 
original  version  of  the  code  (1.10--10/30/84)  did  not  implement  a  physically  correct  initial 
condition,  but  rather  set  the  initial  temperature  of  the  substrate  equal  to  zero.  This 
simplified  approach  was  accomplished  by  setting  all  the  yi's  equal  to  zero  to  start  the 
integration.  The  simplified  initial  condition  caused  a  large  initial  spike  in  the  calculated 
heat-flux  value  which  quickly  decayed  into  the  steady-state  solution.  This  physically 
incorrect  initial  transient  was  corrected  by  implementing  a  more  realistic  initial  condition. 

Non-uniform  initial  conditions  are  considered  because  it  is  often  desirable  to  only 
reduce  a  portion  of  the  recorded  data.  If  the  data  reduction  begins  at  a  time  after  the  tunnel 
flow  begins,  the  substrate  is  no  longer  at  a  uniform  temperature.  A  non-uniform  initial 
condition  can  be  estimated  for  a  finite  thickness  case  because  the  lower  surface  temperature 
is  known.  However,  there  is  no  accurate  method  to  specify  a  non-uniform  initial  condition 
for  the  semi-infinite  thickness  case.  Therefore,  the  following  two  approaches  are  used. 
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2.2.3 .2.1  --  Finite  Thickness  Circuit  Initial  Conditions 

To  account  for  non-uniform  initial  conditions,  a  linear  distribution  between  the 
known  upper  and  lower  surface  temperatures  is  assumed.  The  linear  distribution 
corresponds  physically  to  steady  state  heat  conduction  through  the  gauge  substrate. 

T(x)  =  T„+(T..,-To)i 

where  Tq  is  the  upper  surface  temperature  at  the  start  time  and  Tn+l  is  the  bottom  surface 
temperature  at  the  start  time.  In  terms  of  the  analogue. 


Therefore,  the  initial  voltage  at  any  circuit  node  "m"  is 


t 

V.=V„  +  (V,.,-V„)R,-2^r;  (2.40) 

i=l 

where  m  =  1,  2,  3,  . . .,  n. 

These  initial  voltages  are  then  transformed  into  y  coordinates  to  start  the  integration. 
Recall  the  first  transformation  to  a  synunetric  coefficient  matrix  given  by  Equation  2.14. 

u,  =  v,V^ 

The  second  transformation  to  the  y  coordinate  is  given  by  Equation  2.9. 

y  =  [Qru 

where  [Q]'^  is  the  transformed  matrix  of  eigenvectors.  The  transformation  from  a  voltage 
at  node  "m"  to  the  corresponding  ym  is  given  by 
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(2.41) 


i=l 


where  again  m  =  1,  2,  3, . . n. 


2.2.3.2.2  ~  Semi-Infinite  Thickness  Circuit  Initial  Conditions 

For  the  semi-infinite  case  the  integration  must  be  started  at  a  time  when  the  initial 
temperature  distribution  in  the  substrate  is  accurately  known.  For  all  practical  cases  this  is 
at  a  time  before  the  experiment  is  initiated  when  the  substrate  is  at  a  known  uniform 
temperature.  In  this  case  the  initial  y's  are  given  as 

n  _ 

(2.42) 


2.23.2.3  —  Integration  of  ODE 


All  the  ODE'S  are  integrated  over  each  time  interval  in  the  upper  surface  temperature 
file.  For  the  finite  case,  a  lower  surface  temperature  is  read  at  the  same  time.  The 
integration  is  carried  out  with  a  fourth  order  Runge-Kutta  method  where 

f=f(t.y) 

The  four  steps  are: 


An  =  f(t,y)dt 


„  J  dt  An'\ 
Bn  =  f  t+—,y  +  — 

V  2  2  J 


^  J  dt  Bn^ 
Cn  =  f  X  +—,y  +  — 

\  2  2  J 


and  finally. 


Dn  =  f(t-t-dt,y-i-Cn) 
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This  is  coded  as 


y(t  +  dt)  =  y(t)  + 


An  +  2Bn  +  2Cn  +  Dn 

6 


*_  dt  ,  ,  quVo(t)  ,  qn.iVn.l(t) 

All;  =  - i  Ajy;  + - :  . i"-:—  H - ; - 7=7=- 

R,CJ  “ 


Bn;  = 


dt 


‘  RiQ 


lu 


V„(t)+V„(t+dt)' 

2 


r;V^ 


qn.i 


Vn.l(t)  +  Vn.7(t  +  dt)' 
2 


Rn+lV^n 


Cn;  = 


dt 


■  RiC, 


Bn  A 

yi+-r- 

2  y 


flu 

rv„(t)+v„(t+dt)i 

n 

rv„,(t)+v..,(t+dt)i] 

2 

‘  -f- 

2  J 

R,V^ 

^n+lV^ 

Dn;  = 


dt 


‘  RiQ 


RiVq  r„.iV^  J 


Thus  giving,  y i  (t  +  dt) 

The  final  task  is  to  find  the  voltage  at  node  one  to  evaluate  the  heat  flux.  Working 
back  through  the  transformations, 

n 

u,  =  ^quyi 

i=l 

Vi  = 

Then,  the  surface  heat-flux  value  is  given  by 

V  -V 

q  =  I=  Q— ^  (2.43) 

R. 
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SECTION  3  -  CODE  VERMCATION 


Although  the  data  reduction  code,  known  as  ACQ,  has  been  in  use  for  nearly  ten 
years  at  MIT,  few  details  are  available  on  its  verification.  Therefore,  a  code  verification 
was  conducted  at  Calspan  to  document  the  code's  capabilities  and  accuracy. 

Epstein,  Guenette,  Norton  and  Yuzhang  [1985]  reported  that  the  accuracy  of  ACQ 
was  established  by  back  calculating  the  surface  heat-flux  for  exact  analytical  upper  and 
lower  surface  temperature  solutions  in  two  cases:  a  step  change  in  heat-flux,  and  a 
sinusoidal  heat-flux  variation.  They  used  this  same  approach  to  evaluate  the  number  of 
nodes  required  for  a  specific  bandwidth  and  accuracy.  It  was  found  that  nine  nodes  yielded 
agreement  to  better  than  1%  over  the  entire  DC  to  100  kHz  fi'equency  domain  including  less 
than  0.3%  error  in  the  5-2000  Hz  band. 

The  code  verification  carried  out  here  is  two  fold.  The  first  part  revisits  the  two  test 
cases  studied  by  Epstein  et  al.  [1985]  and  adds  a  third  semi-infinite  case.  All  the  cases  use 
the  recommended  nine  nodes.  The  second  part  of  the  verification  process  illustrates  how 
anomalies  in  the  recorded  signals  affect  the  accuracy  of  the  code's  solution.  The  two 
anomalies  studied  are:  digitization  error  due  to  low  analog  to  digital  resolution,  and  random 
noise  super-imposed  onto  the  signal.  The  verification  process  concludes  with  a 
comparison  of  ACQ  with  Cook/Felderman  [1966]  using  actual  data  from  a  semi-infinite 
Pyrex  gauge  located  on  the  SSME  turbine  stage. 
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3.1  -  Exact  Analytical  Test  Cases 


Epstein  et  al.  [1985]  provide  analytical  upper  and  lower  surface  temperature 
solutions  for  the  case  of  a  step  change  in  surface  heat-flux,  and  for  the  case  of  a  steady 
harmonic  variation  in  surface  heat-flux.  These  exact  temperature  profiles  are  used  as  input 
to  the  code  for  back  calculating  the  surface  heat-flux.  A  third  test  case  taken  from  Incropera 
and  DeWitt  [1985]  for  convective  heat  transfer  at  the  surface  of  a  semi-infinite  slab  is 
performed  for  verification  of  the  semi-infinite  calculations. 


3.1.1  -  Steady  Harmonic  Surface  Heat-flux  Variation 


The  analytical  temperature  distribution  through  a  gauge  substrate  subjected  to  a 
steady  harmonic  surface  heat-flux  is  given  below.  The  surface  heat-flux  is  taken  to  be 

=  Q -1- Qo cos(oot)  for  -«»<t<-(-oo  (3.1) 

x=0 

where  Q  is  the  mean  or  DC  heat-flux  level;  Qo  is  the  perturbation  in  heat-flux  about  the 
mean,  and  co  is  the  frequency  of  the  perturbation.  The  mean  heat-flux  is  given  by 

Q  =  ^(T(0)-T(d))  (3.2) 

where  T(0)  and  T(d)  are  the  mean  upper  and  lower  surface  temperatures  respectively. 

The  temperature  distribution  in  the  substrate  is 

T(x,t)  =  T(x)  +  QoM— cos(©t-(t))  for  0<x<d  (3.3) 

k 


where. 


M  =  V2 


I  ■ 

2 

2  X 

I“cj 

exp 

d 

U^  +  D^J 

CO 

v“cy 


x  % 
—  +  — 


d  4 


-Arc  tan 


BC-ADA 

ac+bdJ 
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A  sensor  characteristic  frequency  (o^  is  introduced  based  upon  the  substrate  thickness  d, 
which  is  defined  as 


“'”2UJ  pCpk 


The  other  terms  in  Equation  3.3  are  defined  as  follows: 


Three  test  cases  are  used  to  study  ACQ's  accuracy  in  calculating  the  steady 
harmonic  surface  heat-flux.  The  influence  of  the  sample  frequency  and  the  influence  of  the 
signal  frequency  are  ascertained  with  the  two  additional  tests.  The  input  parameter  values 
for  the  three  test  cases  are  given  in  Table  3.1.  The  second  test  differs  from  the  first  only  in 
its  sample  frequency  and  the  third  test  differs  from  the  first  only  in  its  blade  passing 
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frequency.  The  values  chosen  are  representative  of  a  typical  Kapton  gauge  used  to  measure 
the  heat-flux  in  a  turbine  tested  within  a  short  duration  facility. 


(pCp>=),„„ 

330,625. 

(W's)/(m'‘K') 

69,334,420. 

(w"s)/(m"K") 

k/d 

8,486 

W/(m'^K) 

Q 

851,250 

W/m^ 

T(0) 

400 

K 

T(d) 

300 

K 

Qo 

681,000 

W/m^ 

(0 

50,265 1/s  (blade  pass  freq.  of  8,000  Hz),  tests  1  &  2 

1/s  (blade  pass  fireq.  of  4,000  Hz),  test  3 

fs 

100 

kHz  tests  1  &  3 

200  kHz  test  2 


25,132 


Figure  3.1  gives  the  analytical  upper  and  lower  surface  temperatures  input  into 
ACQ  for  test  number  1.  Figure  3.2  is  a  plot  of  both  the  theoretical  heat-flux  and  the 
calculated  ACQ  heat-flux  versus  time  for  test  number  1.  This  plot  shows  that  the  ACQ 
heat-flux  lags  the  theoretical  by  approximately  0.005  ms  which  is  one  half  a  sample  period. 
However,  Figure  3.2  also  shows  that  the  calculated  heat-flux  is  nearly  identical  to  the 
theoretical  heat-flux  in  every  other  respect;  i.  e.,  frequency,  perturbation  magnitude  and 


mean. 

The  phase  lag  in  the  ACQ  calculations  was  first  investigated  under  the  assumption 
that  it  was  a  programming  error.  However,  no  error  could  be  identified.  Closer 
examination  suggested  that  it  might  be  a  numerical  inaccuracy.  The  problem  can  be 
visualized  by  considering  the  theoretically  phase  lag  between  the  temperature  and  heat-flux 
given  by  ({)  in  Equation  3.3.  It  appears  that  this  phase  lag  is  not  being  accurately  duplicated 
by  ACQ.  If  the  problem  is  numerical,  the  error  should  be  reduced  by  a  decrease  in  the 
integration  step  size.  This  hypothesis  is  investigated  by  test  number  2. 
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Doubling  the  sample  frequency  in  test  2  reduces  the  integration  step  size  by  one 
half.  The  result  for  this  change  is  given  in  Figure  3.3  as  a  plot  of  the  theoretical  and 
calculated  heat-flux.  The  lag  between  the  calculated  and  theoretical  heat-flux  is 
approximately  half  that  in  test  1.  As  in  test  1,  the  temporal  1^  is  approximately  half  a  time 
step. 

From  a  practical  stand  point  of  reducing  turbine  data,  the  phase  error  is 
inconsequential  as  long  as  all  the  frequencies  in  the  signal  are  shifted  by  the  same  amount. 
To  ensure  that  this  is  the  case,  test  3  was  conducted  which  has  half  the  blade  passing 
frequency  of  test  1.  Figure  3.4  is  a  plot  of  the  theoretical  and  calculated  heat-flux  for  test 
number  3.  These  results  show  the  same  phase  lag  as  test  1  indicating  that  the  lag  is  not  a 
function  of  the  signal  frequency. 

To  determine  the  accuracy  of  the  steady  harmonic  heat-flux  calculations,  the 
theoretical  heat-flux  was  shifted  half  a  time  step  into  phase  with  the  calculated  heat-flux  to 
determine  a  meaningful  difference  between  the  two.  Figure  3.5  gives  the  phase  shifted 
theoretical  heat-flux  and  the  calculated  heat-flux  along  with  the  difference  between  the  two 
for  test  number  1.  The  plot  shows  that  the  calculated  heat-flux  is  larger  than  theory  while 
the  function  is  increasing  and  lower  while  it  is  decreasing.  The  largest  errors  occur  just 
after  the  min.  and  max.  inflections.  Figures  3.6  and  3.7  give  the  absolute  and  percent 
differences  respectively  for  test  number  1.  Several  points  due  exceed  the  1%  NUT  stated 
accuracy  in  this  test  case. 


3.1.2  -  Step  in  Surface  Heat-flux 

In  the  step  change  in  surface  heat-flux  case  both  the  transient  and  steady  state 
solutions  to  the  problem  are  of  interest.  The  initial  condition  is 

T(x,0)  =  T.  (3.4) 
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For  the  surface  heat-flux, 


.  ,31  fO,-°o<t<to 

=  “k-— -  =  ^ 

5xLo  lQo.to^t<~ 

the  temperature  distribution  throu^  the  substrate  is  given  by 


(3.5) 


where. 


T(x.t)  =  T,  +Mrij|^(_r)“[ierfc(Y.)-nerfc(Z.)] 

^  n=0 

x/d  2n 


(3.6) 


_  x/d  2(n-i-l) 

and  ierfc( )  is  the  first  integral  of  the  complementary  error  function.  A  gauge  time  constant 
X  is  introduced  based  upon  the  gauge  thermal  diffiisivity  and  substrate  thickness. 


X 


(3.7) 


The  input  parameter  values  for  this  test  case  are  given  in  Table  3.2.  Again  they  are 
representative  of  a  typical  Kapton  gauge. 


Table  3.2  —  Input  Parameters.  Step  Change  in  Surface  Heat-flux 


(pCpk) 

V  P  /gauge 

330,625. 

(W's)/(m^K^) 

fpCpk) 

p  /blade 

69,334,420. 

(W's)/(m''K-) 

k/d 

8,486 

W/(m-K) 

Qo 

851,250 

W/m- 

to 

5.0 

ms 

Ti 

300 

K 

fs 

100 

kHz 
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The  analytical  upper  and  lower  surface  temperatures,  given  by  Equation  3.6,  are 
shown  in  Figure  3.8  for  the  parameters  of  Table  3.2.  These  temperatures  produce  the 
calculated  step  in  heat-flux  shown  in  Figure  3.9.  Figure  3.10  better  shows  the  initial  over 
shoot  in  the  calculated  heat-flux.  The  absolute  and  percent  differences  between  the 
calculated  and  theoretical  heat-flux  steps  are  given  in  Figures  3.11  and  3.12  respectively. 

These  results  show  that  the  calculations  recover  in  less  than  0.3  ms  (30  time  steps) 
to  match  the  magnitude  of  the  step  almost  exactly.  The  over  shoot  in  the  calculated  heat- 
flux  is  a  numerical  error  due  to  the  upper  surface  temperature's  abrupt  change  in  slope  and 
the  very  large  initial  value  of  that  slope.  The  4%  error  of  the  over  shoot  is  undesirable  but 
acceptable  considering  its  short  duration  along  with  the  fact  that  a  step  change  is  a  worse 
case  example. 


3.1.3  “  Convective  Heat  Transfer  at  the  Surface  of  a  Semi-infinite  Slab 


This  test  case  examines  the  heat  transfer  at  the  surface  of  a  semi-infinite  slab  given 
the  initial  condition 


T(x,0)  =  Ti 

(3.8) 

and  boundary  conditions 

Interior: 

T(oo,t)  =  T; 

(3.9) 

Surface: 

-k^ 

ax 

=  h[T„-T(0,t)] 

x=0 

(3.10) 

where  h  is  the  convective  heat  transfer  coefficient  and  T„  is  the  free  stream  fluid 
temperature.  The  closed  form  solution  for  the  temperature  distribution  in  the  slab  as  given 
by  Incropera  and  DeWitt  [1985] 


T(x,t)-T, 

T„-T. 


f 


=  erfc 

V 


X 

2'\/t 


V 


h^t 

pCpky 


r 


erfc 


X 

iVT 


(3.11) 
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The  input  parameter  values  for  this  test  case  are  given  in  Table  3.3.  These  values 
are  representative  of  a  Pyrex-type  gauge. 

Table  3.3  -  Input  Parameters.  Semi-infinite  Slab  (Pvrex  Insert  Gauge) 


2.84x106 

(w's)/(m"K") 

k/d 

551 

W/(m^K) 

d 

0.00254 

m 

Ti 

294 

K 

T„ 

533 

K 

h 

1260 

W/(m"K) 

fs 

20 

kHz 

The  analytical  upper  surface  temperature  generated  by  Equation  3.11  is  given  in 
Figure  3.13.  The  heat-flux  calculated  with  this  temperature  trace  is  given  along  with  the 
theoretical  heat-flux  in  Figure  3.14.  The  absolute  and  percent  differences  between  the 
calculated  and  theoretical  heat-flux  are  given  in  Figures  3.15  and  3.16  respectively.  These 
results  show  an  initial  error  of  nearly  10%  in  the  calculations  which  disappears  in 
approximately  0.3  ms  to  settle  at  a  value  near  0.1%. 

3.2  —  Effects  of  Signal  Anomalies 

It  is  necessary  to  understand  how  errors  in  the  recorded  signals  propagate  through 
the  heat-flux  calculations.  Two  typical  sources  of  signal  errors  are  studied  here.  The  first 
source  considered  is  analog  to  digital  conversion  (A/D)  which  introduces  error  because  of 
the  finite  resolution  of  A/D  systems.  The  second  source  is  random  noise  which  is  super¬ 
imposed  onto  the  desired  signal. 
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3.2.1  —  A/D  Resolution  Error  Propagation 


The  effects  of  finite  A/D  resolution  are  studied  using  the  steady  harmonic  surface 
heat-flux  case  with  the  same  parameters  as  in  test  number  1.  The  study  simulates  the  DC 
coupled  acquisition  of  the  temperature  signals  which  have  small  perturbations  in 
comparison  to  their  mean  levels;  reference  Figure  3.1.  Therefore,  the  range  of  perturbation 
temperatures  is  characterized  by  a  small  percentage  of  the  A/D  system's  available  bits.  In 
the  case  presented  here,  1.5%  of  4096  available  bits  were  used  to  characterize  the  +/-  3  K 
upper  surface  temperature  perturbation  about  its  400  K  mean  value. 

The  calculated  and  theoretical  heat-flux  plots  for  this  test  case  are  shown  in  Rgure 
3.17.  This  figure  appears  to  be  identical  to  Figure  3.2,  but  the  differences  given  in  Figures 
3.18  and  3.19  show  the  truncation  errors  caused  by  the  low  A/D  resolution.  The  error  for 
this  case  is  as  high  as  5%  of  the  mean  heat-flux  and  13.5%  of  the  local  heat-flux.  This 
represents  a  factor  of  three  increase  in  the  error  strictly  due  to  the  low  A/D  resolution. 

The  test  demonstrated  here  is  a  worse  case.  This  type  of  A/D  truncation  error  can 
be  gready  reduced  by  utilizing  the  A/D  system's  resolution  to  its  fullest.  For  instance,  a 
DC  offset  can  be  used  to  zero  the  smallest  expected  voltage  reading.  Then,  a  gain  can  be 
applied  to  the  expected  signal  so  that  it  will  fill  the  full  A/D  range.  The  signal  should  be 
acquired  in  a  DC  couple  mode  because  of  the  decaying  nature  of  the  signal  in  short  duration 
test  facilities.  However,  AC  coupled  acquisition  can  be  performed  over  short  (quasi 
steady)  periods  with  the  DSP  data  channels  but  is  not  recommended  with  the  Data  Lab 
channels  because  of  their  poor  RC  time  constant. 

3.2.2  —  Random  Noise  Error  Propagation 

The  steady  harmonic  surface  heat-flux  test  case  is  used  to  study  how  random  signal 
noise  affects  the  accuracy  of  the  heat-flux  calculations.  Signal  noise  is  simulated  by  a 
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random  number  generator  which  super-imposes  an  error  onto  the  exact  temperature  signal. 
This  random  error  is  applied  to  the  signal  at  the  digital  level  in  terms  of  plus  or  minus  a 
number  of  bits.  The  noise  has  a  maximum  possible  value  which  is  specified  as  a 
percentage  of  the  full  A/D  bit  range. 

For  the  case  studied  here,  the  A/D  system  is  assumed  to  have  a  4096  bit  range.  The 
exact  temperatures  from  test  number  1  are  characterized  with  1.5%  of  the  available  bits  in 
the  simulated  DC  couple  mode.  A  0.2%  of  full  scale  noise  level  is  super-imposed  onto  the 
temperatures.  This  level  of  noise  corresponds  to  approximately  +/-  0.3  K  for  this  particular 
test  case.  Figure  3.20  compares  the  exact  upper  temperature  from  Figure  3.1  to  the 
simulated  upper  temperature  with  the  A/D  resolution  error  and  random  noise.  The  signals 
only  differ  slightly  with  the  effect  of  the  random  noise  being  seen  in  the  separation  of  some 
of  the  data  points. 

The  ACQ  heat-flux  calculated  from  the  simulated  temperature  data  is  given  in  Figure 
3.21  along  with  the  theoretical  heat-flux.  Again  the  phase  shift  is  present  in  the  calculations 
and  the  effect  of  the  signal  anomalies  is  discernible.  Figures  3.22  and  3.23  give  the 
absolute  and  percent  differences  which  indicate  the  random  nature  of  the  error.  While  most 
of  the  calculations  are  within  +/-  5%,  there  are  a  few  that  are  over  20%. 

This  test,  like  the  previous  one,  is  a  worse  case  scenario.  For  example,  the  A/D 
resolution  is  much  lower  than  could  be  achieved,  as  explained  in  Section  3.2.1.  Also, 
Epstein  et  al.  [1985]  states  that  electrical  signal-to-noise  considerations  are  not  important 
and  quotes  an  analysis  performed  by  Guenette  [1985]  which  shows  that  the  noise 
equivalent  temperature  is  only  0.002  K. 
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3.3  —  Actual  Data  Reduction 


The  last  part  of  the  code  verification  process  entails  the  reduction  of  heat-flux  gauge 
data  from  a  turbine  experiment.  Data  firom  a  single  sensor  semi-infinite  Pyrex  gauge  was 
chosen  because  it  could  be  accurately  reduced  by  an  uiuelated  method  developed  by  Cook 
and  Felderman  [1966].  The  Cook/Felderman  method  is  regarded  as  a  tried  and  true  method 
for  reducing  semi-infinite  heat-flux  gauge  data.  The  heat-flux  data  was  taken  firom  an 
SSME  turbo  pump  test  conducted  in  Calspan's  Turbine  Test  Facility  (TTF).  The  recorded 
gauge  surface  temperature  is  given  in  Figure  3.24. 

The  heat-flux  calculated  by  both  ACQ  and  Cook/Felderman  are  given  m  Rgure 
3.25.  The  plot  shows  that  the  two  methods  are  in  excellent  agreement.  A  more  detailed 
plot  over  the  usable  test  time  is  given  by  Figure  3.26.  This  plot  shows  that  the  ACQ 
calculations  are  consistently  2.0%  to  2.5%  lower  than  the  Cook/Felderman  calculations. 
Considering  the  data  and  the  difference  in  the  methods,  the  correlation  is  very  encouraging. 

3.4  —  Results  Summary 

•  Steady  harmonic  surface  heat-flux:  +3.25%  to  -1.2%  error 

•  Step  change  in  surface  heat-flux:  +3.25%  (initial  over  shoot  for  -0.3  ms)  to  +0.05% 

•  Convective  heat  transfer  over  semi-infinite  slab:  +10%  (initial  error  for  -0.3  ms)  to 
+0.1%  error 

•  A/D  resolution  error:  +13.5%  (most  points  below  +5%)  to-3.2%  error 

•  A/D  resolution  plus  0.2%  full  scale  random  noise  error:  +28%  to  -22%  error  (most 
points  between  +/-  5%) 

•  Pyrex-type  gauge  data  from  SSME  turbine  test:  -2.5%  to  -2.0%  difference  between 
Cook/Felderman 
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3.5  —  Conclusions 


The  ACQ  heat-flux  data  reduction  code  has  been  verified.  The  code  provides 
physically  correct  and  accurate  results  when  compared  to  various  exact  analytical  solutions 
and  one  other  accepted  data  reduction  code.  The  propagation  of  signal  error  through  the 
data  reduction  calculations  has  been  investigated.  Under  what  should  be  worse  case  data 
acquisition  conditions,  the  simulated  data  errors  produced  significant  but  acceptable  errors 
in  the  heat-flux  calculations.  The  random  errors  encountered  would  be  nearly  eliminate 
with  adequate  ensemble  averaging. 
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Figure  3.1  --  Exact  analytical  upper  and  lower  surface  temperatures  for  steady 
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Figure  3.4  —  Theoretical  and  ACQ  calcul 
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Figure  3.6  ~  Absolute  difference  between  theoretical  and  ACQ  steady  harmonic  surface  heat-flux 
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Figure  3.8  -  Exact  analytical  upper  and  lower  surface  temperatui 
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Figure  3.10  --  Theoretical  and  ACQ  calculated  step  change  in  surface 
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Figure  3.1 1  —  Absolute  difference  between  theoretical  and  ACQ  step  change  in  surface  heat-flux 
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Figure  3.12  -  Percent  difference  between 
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Figure  3.14  —  Theoretical  and  ACQ  calculated  heat-flux  for  a  semi-infinite  slab  in  convection 
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Figure  3.17  --  Theoretical  and  ACQ  calculated  steady  harmonic  heat-flux  with 
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Figure  3.18  --  Absolute  difference  between  theoretical  and  ACQ  steady  harmonic  heat-flux  with  1.5%  of  the  A/D  resolution 
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Figure  3.21  -  Theoretical  and  ACQ  calculated  steady  harmonic  surface  heat-flux  with  1.5%  A/D  resolution  and  0.2%  rai 
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Figure  3.22  —  Absolute  difference  between  theoretical  and  ACQ  steady  harmonic  surface  heat-flux  with  1.5%  A/D  resolution  and 

0.2%  random  noise 


Figure  3.23  --  Percent  difference  between  theoretical  and  ACQ  steady  harmonic  surface  heat-flux  with  1.5%  A/D  resolution  and 

0.2%  random  noise 
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Figure  3.26  -  ACQ  and  Cook/Felderman  calculated  heat-flux  for  Pyrex-type  gauge  on  SSME  turbine 


SECTION  4  -  CODE  APPLICATION 


This  section  provides  the  user  with  the  information  necessary  to  properly  apply 
ACQ.  ACQ  is  compatible  with  the  Data  Acquisition  System  (DAS).  It  operates  with  DAS 
format  input  and  output  files  reading  and  creating  the  appropriate  file  headers.  The 
program  can  be  operated  in  two  modes:  interactive  and  automatic.  The  interactive  mode 
provides  the  user  with  a  convenient  tool  for  generating  heat-flux  time  histories  for  one  or 
two  gauges.  The  automatic  mode  is  more  like  a  batch  process  which  not  only  performs 
multiple  gauge  heat-flux  calculations,  but  also  places  pertinent  information  into  the  mn  log 
files. 

The  first  part  of  this  section  outlines  the  code's  two  modes  of  operation  with 
emphasis  on  the  required  input.  The  second  section  examines  some  important 
considerations  in  using  ACQ.  These  considerations  point  out  that  the  data  reduction 
process  must  be  considered  in  the  design  of  the  experiment. 

4.1  -  Operation  of  ACQ 

ACQ  is  fully  automated  so  that  multiple  gauge  heat-flux  calculations  can  be  invoked 
with  a  single  command  line.  In  this  automatic  mode,  the  code  operation  is  controlled 
through  a  sensor  details  file.  As  indicated  earlier,  the  automatic  mode  can  be  passed  by  to 
give  the  user  greater  versatility  in  an  interactive  mode.  However,  the  sensor  details  file  can 
still  be  utilized  within  the  interactive  mode. 
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4.1.1  --  Sensor  Details  File 


The  sensor  details  file  is  a  user  created  listing  of  similar  gauges  that  are  all  mounted 
on  a  given  test  article.  For  instance,  a  sensor  details  file  would  typically  contain  a  listing  of 
all  the  double-sided  gauges  on  a  particular  turbine  with  a  line  in  the  file  for  each  test  run 
recorded.  Each  line  contains  information  on  the  gauge  and  the  run  data.  The  user  is  free  to 
break  up  the  gauges  and  runs  into  as  many  sensor  files  as  is  convenient. 

A  line  in  the  sensor  details  file  contains  the  gauge/run  label  in  single  quotes,  the 
gauge  span  and  wetted  distance  coordinates,  the  pCpk  [W^s/m^k]  and  the  k/d  [W/m^k]  of 
the  gauge.  The  gauge/run  label  is  the  six  character  DAS  designation  for  the  particular  test 
run  and  sensor.  This  label  is  used  as  a  prefix  for  the  input  temperature  file  names;  i.e., 
'LABEL'T.l  for  top  temperature  file  and  'LABEL’S.  1  for  bottom  temperature  file. 
Likewise,  the  heat-flux  output  file  generated  by  ACQ  is  'LABEL'Q.2. 

4.1.2  —  Interactive  Mode 

Operation  of  the  ACQ  interactive  mode  is  outlined  below  with  an  explanation  the 
various  prompts. 

•  Type  'acq'  to  enter  interactive  mode. 

•  'PROCESS  A  SERIES  OF  DATA  FILES?  (Y/N)'. 

'Y':  Enter  the  sensor  details  file  name.  Program  then  allows  the  user  to  delete  any 
items  in  the  file  not  wanted  for  processing. 

'N':  Enter  gauge  label,  pCpk,andk/d.  Program  then  allows  the  user  to  edit  any 
of  the  three  inputs  before  beginning. 

•  'HNITE  LENGTH  (0)  OR  SEMI-INHNITE  ( 1 )  LINE  (0/1 )  ?' 

'O'  for  double-sided  gauge 
'!'  for  single-sided  gauge 
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•  If  only  one  gauge  is  being  processed,  the  program  provides  the  versatility  to  use 
filenames  other  than  the  displayed  defaults. 

'ENTER  FILENAME  FOR  FIRST  UPPER  TEMP.  [DEF  'LABEL’T.l]' 

'ENTER  FILENAME  FOR  FIRST  LOWER  TEMP.  PEF  'LABEL'S.  1]’ 

'ENTER  FILENAME  FOR  FIRST  HEAT  FLUX  PEF  'LABEL'Q.2]' 

•  'ENTER  MAXIMUM  FREQUENCY  (HZ)  FOR  ANALOGUE  PEF=fs/2]'  The 
default  (1/2  sample  frequency)  is  the  largest  possible  frequency  allowable  before 
frequency  folding  occurs.  A  frequency  less  than  the  default  may  be  entered. 

•  'ENTER  REQUIRED  NUMBER  OF  STAGES  [DEF=9]'  Experience  indicates  that  9 
analogue  nodes  is  optimum  for  the  program.  However,  for  unusual  circumstances  this 
may  be  changed. 

•  'CAPACITOR  GEOMETRIC  VALUE  IS  '  y'.  IS  THIS  ACCEPTABLE  (Y/N)? 
PEF=Y]'  ACQ  allows  the  user  to  review  the  calculated  geometric  factor.  If  for  some 
reason  the  value  is  deemed  unacceptable,  the  user  can  change  the  number  of  stages  until 
the  geometric  factor  is  acceptable. 

•  If  a  finite  length  line  is  used,  the  program  prompts  'ENTER  START  TIME  (MS)'.  The 
start  time  for  a  semi-infinite  line  is  automatically  zero. 

•  'ENTER  END  TIME  (MS)'.  This  can  be  less  than  or  equal  to  the  end  of  the  data  file. 

•  'SCREEN  OUTPUT  EVERY  N  STEPS,  ENTER  N'.  ACQ  allows  the  user  to  view 
the  calculated  heat-flux  values  on  the  screen. 

•  ACQ  then  creates  a  DAS  output  file  of  the  specified  name  with  the  appropriate  header 
information  and  heat-flux  values. 


4.1.3  —  Automatic  Mode 

The  automatic  mode  allows  ACQ  to  be  operated  from  the  Data  Reduction  Process 
(DRP)  environment  through  an  'XU'  interface.  ACQ  is  invoked  with  the  following 
command  line:  'acq  FILENAME  [STARTTIMEpNDTIME[I_SINF]]]'. 

•  'FILENAME'  refers  to  the  sensor  details  file  that  contains  the  gauge/ran  labels  to  be 
reduced. 

•  'STARTTIME'  is  in  ms.  For  a  semi-infinite  analogue,  this  is  set  to  zero,  and  if  omitted 
the  default  is  zero. 

•  'ENDTIME'  is  in  ms.  If  omitted,  the  default  is  the  end  of  the  data  file. 
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•  'I_SINF  is  an  integer  indicator  to  select  the  semi-infinite  analogue.  If  omitted  or 
zero,  the  finite  analogue  is  used;  otherwise,  the  semi-infinite  analogue  is  used. 

When  ACQ  is  accessed  from  the  DRP,  the  user  is  informed  from  the  following  list 
of  error  messages  of  any  error  which  causes  the  code  to  abort; 

1. )  sensor  details  file  does  not  exist 

2. )  error  reading  sensor  details  file 

3. )  required  ’INDEXER.!’  file  does  not  exist 

4. )  no  data  in  the  specified  time  interval 

5. )  upper  and  lower  temperature  files  recorded  differently 

6. )  start  time  error 

7. )  end  time  error 

8. )  convergence  error  in  eigenvalue  problem 


4.2  —  Data  Reduction  Considerations 

This  section  reemphasizes  some  of  the  subtle,  yet  important,  considerations  for 
applying  ACQ.  For  instance,  for  a  finite  analogue  the  upper  and  lower  surface 
temperatures  must  correspond  to  the  same  time;  i.e.,  simultaneously  sampled.  Also  the 
finite  analogue  allows  the  user  to  reduce  the  data  over  a  smaller  subinterval  of  the  total  test 
time  by  choosing  a  start  time  other  than  zero. 

In  applying  ACQ  with  a  semi-infinite  analogue,  there  are  two  important 
considerations  to  review.  The  first  is  the  semi-infinite  assumption  which  must  be  valid  for 
the  data  being  reduced.  This  assumption  is  met  if  the  test  time  is  less  than  the  reciprocal  of 
the  minimum  frequency  calculated  according  to  Equation  2.30;  i.e.. 


fend  f  start 


< 


1 


UJ  pCk 


(4.1) 
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The  second  consideration  with  a  semi-infinite  analogue  is  the  initial  condition. 

As  stated  in  Section  2. 2.3. 2.2,  the  initial  condition  must  be  a  known  uniform  substrate 
temperature.  This  is  accomplished  by  ensuring  that  the  first  data  point  (tstart  =  0)  is  prior 
to  any  test  section  flow.  The  next  very  important  consideration,  discussed  in  more  detail 
below,  is  the  stability  of  the  integration  which  applies  to  both  finite  and  semi-infinite 
analogues. 


4.3.1  -  Stability 


The  fourth  order  Runge-Kutta  integration  method  used  by  ACQ  is  conditionally 
stable.  Hoffman  [1990]  notes  that  the  following  condition  must  be  met  for  this  method  to 
be  stable. 


a(At)<  2.785  (4.2) 

where  a  is  taken  from  the  homogeneous  Ordinary  Differential  Equation  (ODE)  of  the  form 

^  =  -ay(t) 
dt 

The  ODE  integrated  within  ACQ  is  given  by  Equation  2.13.  Its  corresponding 
homogeneous  ODE  is 


thus. 


a  =  ■ 


The  stability  criterion  for  ACQ  is  therefore  given  by 

-X 


f  >■ 


(2.785)RiC, 

where  X  is  the  largest  negative  eigenvalue  for  the  analogue's  system  of  equations. 


(4.3) 
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Examination  of  Equation  4.3  indicates  that  stability  is  a  function  of  the 
following  input  parameters:  pCpk,  k/d,  no.  of  nodes,  and  maximum  frequency. 
Unfortunately,  Equation  4.3  is  not  simple  to  evaluate  because  it  requites  the  eigenvalues  of 
the  analogue's  system  of  equations.  To  illustrate  some  typical  minimum  sample 
frequencies.  Equation  4.3  is  evaluated  for  the  two  gauges  used  in  the  Section  3  verification 
tests. 

The  double-sided  Kapton  gauge  minimum  sample  frequency  for  stability  is 
evaluated  using  the  steady  harmonic  heat-flux  test  number  1.  The  parameters  for  this  test 
are  given  in  Table  3.1,  of  Section  3,  with  fmax  =  50  kHz  and  n  =  9.  The  nine  eigenvalues 
range  from  -1.6418  to  -0.014017.  Therefore,  the  minimum  sample  frequency  for  this  case 
is  89  kHz. 

The  minimum  sample  frequency  for  the  single-sided  Pyrex  gauge  is  evaluated  using 
the  convective  heat  transfer  test  case.  The  parameters  for  this  test  case  are  given  in  Table 
3.3  with  fmax  =  10  kHz  and  n  =  9.  The  eigenvalues  range  from  -1.3548  to  -9.8002xl0‘6. 
Evaluating  Equation  4.3  gives  a  minimum  frequency  of  13  kHz. 

In  both  cases  presented  here,  the  minimum  sample  frequency  for  stabihty  is  less 
than  the  minimum  sample  frequency  required  to  prevent  frequency  folding.  However, 
caution  should  be  taken  in  extending  these  results  to  dissimilar  gauges.  In  fact,  it  might  be 
necessary  to  generate  a  test  case  for  the  particular  gauge  and  sample  frequency  to  ensure 
that  the  code  will  be  stable. 
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program  acq 

^  logarithmic  high-bandwidth  resistance-capacitance 
(R-C)  analogue  to  obtain  high  frequency  heat  transfer  rates. 

Theory  from: 

Oxford  University  Engineering  Laboratory  Report  1382/81 
Authors:  M. L . G. Oldfield,  H.J.Burd  and  N.G.Doe 
"Design  of  Wide-Bandwidth  Analogue  Circuits  for  Heat  Transfer 
Instrumentation  in  Transient  Tunnels" 

Version  1.10  R.J.G. Norton . Oct-30-1984 

Modified  so  that  user  specifies  number  of  stages  and  the 
capacitor  geometric  factor  is  calculated. 

Further  modified  (Fall  1993)  to  correct  initial  boundary 
conditions  and  enable  I/O  of  DAS  format  files. 

The  program  has  been  modified  to  operate  in  two  modes: 
interactive  and  automatic.  The  interactive  mode  provides  the 
user  with  a  convenient  tool  for  generating  heat  flux  time 
histones  for  one  or  two  sensors.  It  is  run  simply  by  enterina 
the  co^and  "acq"  and  supplying  requested  information. 

The  automatic  mode  is  more  like  a  batch  process  and  in 
addition  to  the  heat  flux  attends  to  placing  pertinent  information 

turbine  run  log  files.  The  form  of  the  command  to  execute  the 
program  in  automatic  mode  is:  une 

"acq  filename  [  mintim  [  maxtim  [  i  sinf  ]  ]  ]  " 

Where,  filename  refers  to  the  required  lensor  details  file;  mintim 
and  maxtim  are  respective  min  &  max  times  in  msecs;  and,  i  sinf  is 
an  integer  indicator  to  select  the  semi-infinite  length  analogue 
^^ssing  or  zero,  the  finite  length  analogue  is  used! 
Otherwise,  the  semi-infinite  assumption  is  used  and  mintim  =  zero 
The  sensor  details  file  consists  of  one  or  more  lines  containing, ' 

Single  quotes;  x,  y  coordinates;  and  rho*c*k  ^ 
(W-2-s/mM-K)  &  k/d  (W/m-2-K)  of  the  substrate  material  ?he 

properties  of  the  first  gauge  in  the  sensor  details  files  are  used 
to  determine  global  treatment  of  the  rest.  If  mintim  &  maxtim  are 
entire  recorded  time  interval  of  the  first  gauge  is 

Since  the  data  reduction  process  (drp)  accesses  acq,  there 
must  be  some  way  of  informing  the  user  of  error  conditions.  A 
ist  of  the  conditions  which  cause  the  code  to  abort  follows: 

1.  sensor  details  file  does  not  exist; 

2.  error  reading  sensor  details  file; 

3.,  required_  "INDEXER.  1"  file  does  not  exist; 

4.  no  data  in  the  specified  time  interval; 

5.  upper  and  lower  temp  files  recorded  differently; 

D.  mintim  error; 

7.  maxtim  error; 

8.  minfrq  error; 

9.  convergence  error  in  valvar. 
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I 

Nov 


Work  Arrays  (Size) : 


CAP (MAXST) 

RES (MAXST+1) 

D (MAXST) 

E (MAXST) 

A (MAXST, MAXST) 
Y (MAXST) 

XKOD (MAXPTS) 
RCK (MAXPTS ) 


Capacitance  values  for  each  part  of  the  network 
Resistance  values  for  each  part  of  the  network 
B  matrix  storage  —  diagonal  elements 
B  matrix  storage  —  off— diagonal  elements 
Eigen  vector  storage 

Solution  at  each  node  (transformed  voltaae) 

K/D  for  each  gauge 
RHO*C*K  for  each  gauge 


finclude  <acqchg.inc> 
tinclude  <acqfil.inc> 
tinclude  <acqrpa.inc> 
linclude  <acqsiz.inc> 
tinclude  <dstruc.inc> 

record  /rundes/  trdr,  brdr,  qrdr' 
record  /chdes  /  tcdr,  bcdr,  qcdr 


* 

* 

* 

* 

* 

* 


parameter 

character 

character 


+ 

+ 


integer 

logical 

pointer 

real*4 

real*8 


1 

2 

3 

4 

5 

6 


(maxst  =  40) 

fname*9,  fpos*70,  foutt*70,  foutb*70,  foutq*70 

tin.str'8, 

chirun*6,  ch^out (4 ) *10,  finput*70,  logfil*13 
getlc,  large,  nc_out(4) 
exists,  screen 

(ipvalo,vto) ,  (iptime,  rt) 
vt(*),  vtb(*),  vto(*),  rt(*),  cvalue(4) 
a(40,40),  cap(40),  d(40),  e(40),  res(41),  t(40) 
y(40),  pi,  r8time(128),  sumres,  twopi,  ' 

amp,  an,  beta,  bn,  cl,  cn,  ctot,  deltat,  dgam,  dn,  dt. 

ffunc,  fmax,  fmin,  gam,  gamn,  rl,  rrlcl,  rtot,  voll 
time,  ' 

ul,  vO,  vOb,  vl,  vlb,  v5,  v5b,  vforSb,  vforcl,  vforc5 
vforce,  vforlb,  vforeb  orco, 


Logical  Unit  Numbers: 

1  -  DAS  upper  temperature  file 

2  -  DAS  lower  temperature  file 

3  -  Q  vs  time  output  file 

4  —  Run  log  file 

19  -  Internal  "scratch"  file 


pi  =  dacos (-IdO) 

twopi  =  2d0  *  pi 

call  getenv (' USER' , usrnam) 

open  (19, status=' scratch' ) 

ipvalt  =  0 

ipvalb  =  0 

ipvalo  =  0 

iptime  =  0 

n_args  =  iarge ( ) 

screen  =  .true. 

if (n_args.gt.O)  then 
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screen  =  .false, 
call  getarg (1, f input) 
endif 

Get  descriptor  names,  RHO*C*K  and  K/D 
call  readna ( 0 , 0 , 1 , 1 , npos , screen, f input ) 

Determine  whether  finite  or  semi-infinite  length  analogue 

if  (screen)  then 
isemi  =  -1 
numerr  =  0 

do  while  (isemi.lt. 0  .or.  isemi. gt.l) 

100  write (*,  '  (/a, $) ' ) 

+  'FINITE  LENGTH  (0)  OR  SEMI-INFINITE  (1)  LINE  (0/1)'’  ' 
read(*, *,err=100, end=200)  isemi 
enddo 
endif 

Set  default  data  file  names 

fname  =  name(l) 
ncstem  =  Inblnk (fname) 
foutt  =  fnaime  (: ncstem) //' T.l' 
foutb  =  fname (: ncstem) //'B.l' 
foutq  =  fname (: ncstem) //'Q. 2' 
ncf  =  ncstem  +  3 

Prompt  for  filenames  if  only  one  gauge  is  being  treated 


if  (npos . eq. 1)  then 
numerr  =  0 
if  (screen)  then 

write (*,101)  'UPPER  TEMP. ',  foutt ( :ncf) 
format (/'ENTER  FILENAME  FOR  FIRST  '  ,A,  '  [DEF  '  A  '1') 
read(*, '  (a) '  ,end=200)  fpos  '  ' 

if(fpos.ne.'  ')  foutt  =  fpos 
if (isemi. eq.O)  then 

write  (*,101)  'LOWER  TEMP. ',  foutb ( :ncf) 
read(*, '  (a)' ,end=200)  fpos 
if(fpos.ne.'  ')  foutb  =  fpos 
endif 

write (*,101)  'HEAT  FLUX' , foutq (: ncf ) 
read(*, ' (a) ' ,end=200)  fpos 
if(fpos.ne.'  ')  foutq  =  fpos 
endif 
endif 


Set  up  network  requirements  ********************* 


First:  set  number  of  stages  and  capacitance/resistance  variation 
along  the  line. 


open (1, file  foutt, status—' old' ,  form=' unformatted' 
+  recl=1024,access=' direct') 
read(l, rec=l)  trdr 
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readd,  rec=2)  tcdr 
close (1, status=' keep' ) 

if (screen)  then 
numerr  =  0 

Input  maximum  analogue  frequency 
110  fmax  =  tcdr.lsr  /  2 

write(*,' (a,lpell.4,a,$)') 

+  'ENTER  MAXIMUM  FREQUENCY  (HZ)  FOR  ANALOGUE  rDEF=' 
read(*, ' (a) ',end=200)  fpos 
if(fpos.ne.'  '  .and.  fpos.ne. ' /' )  then 
rewind  19 

(i9f  '  (a) '  )  fpos  ( :  Inblnk  (fpos) ) 
rewind  19 

read(19,*,err=110)  fmax 
if (fmax. le. 0.0)  go  to  110 
endif 
else 

inprun  =  trdr.run 
call  n\im2ch  (inprun,  chirun,nchrun) 
logfil  =  'run'//chirun(:nchrun)//' .log' 
inquire (file=logfil,  exist=exists) 
if (exists)  then 

open (4, file=logfil,  access='  append' ) 
else 

open (4, file=logfil, status=' new' ) 
endif 

call  ufdate (datstr,  timstr) 
fmax  =  tcdr.lsr  /  2 
deltat  =  IdO  /  tcdr.lsr 
St  =0.0 

=  (tcdr.nsampl  -  1)  *  deltat 
isemi  =  0 

if  (n_args.gt. 1)  then 

Second  argument  is  start  time  in  msecs 
call  getarg (2, fpos) 
nfp  =  Inblnk (fpos) 

rewind  19 

write (19,' (a) ' )  fpos (: nfp) 
rewind  19 

fnp0(3:3)  =  char (nfp+48) 
numerr  =  6 

read(19,  fnpO,  err=200)  st 

St  =  St  /  ld3 

if (n_args.gt.2)  then 

Third  argument  is  final  time  in  msecs 
call  getarg (3, fpos) 

=  Inblnk (fpos) 
rewind  19 

write(19, ' (a) ' )  fpos(;nfp) 
rewind  19 

fnp0(3:3)  =  char (nfp+48) 
numerr  =  7 

read(19, fnpO, err=200)  ft 

ft  =  ft  /  ld3 
if (n_args.gt. 3)  then 


, fmax, ' ]  ' 
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Fourth  argument  is  minimum  frequency  in  HZ 
call  getarg ( 4 , fpos ) 
nfp  =  Inblnk(fpos) 
rewind  19 

write  (19, '  (a) ' )  fpos ( :nfp) 
rewind  19 

iw(3:3)  =  char(nfp+48) 
numerr  =  8 

read(19,iw,err=200)  isemi 
if (isemi.ne. 0)  then 
isemi  =  1 

fmin  =  xkod(l)**2  /  (rck(l)  *  twopi) 

St  =  OdO 

endif 
endif 
endif 
endif 

write (4,' (a8, 2x, a9, 2x, alO, 2x, a) ' )  timstr, datstr,usrnam, 

.  'Program  acq  specifics 

ncfin  =  Inblnk (f input) 

write  (4, '  (a) ' )  '*  Sensor  details  taken  from  ' //f input (: ncfin) 
write  (4,'  (a, Ipell . 4, a) ' )  '*  with  RHO*C*K  of  the  substrate  =' , 
+  rck(l),'  W^2-s/m''4-K''2. ' 

if (isemi. eq. 0)  then 

write(4,'  (a, Ipell. 4, a) ' ) 

+  '*  The  finite  length  analogue  with  K/D  =',xkod(l) 

+  '  W/m''2-K  was  used.'  ' 

else 

write (4, ' (a) ' ) 

+  '*  The  semi-infinite  length  analogue  with  K/D  =',xkod(l) 

+  '  W/m''2-K  was  used. '  ' 

endif 

cvalue(l)  =  St  *  1000 
cvalue(2)  =  ft  *  1000 
cvalue(3)  =  fmin 
cvalue ( 4 )  =  f max 
do  i=l,4 

write (fpos, ' (flO . 3) ' )  cvalue (i) 
if (cvalue  (i) . eq.  0 . 0)  then 
ch_out ( i )  =  ' 0 ' 
nc_out ( i )  =1 
else 

11  =  Idblnk(fpos) 

12  =  Inblnk (fpos) 

if (cvalue (i) . eq. aint (cvalue (i) ) )  then 
i2  =  i2  -  4 
else 

do  while  (fpos (i2 : i2 ) . eq. '  0' ) 
i2  =  i2  -  1 
enddo 
endif 

ch_out(i)  =  fpos(il:i2) 
nc_out(i)  =  i2  -  il  +  1 
endif 
enddo 

write  (4, '  (a) ' )  '*  Time  interval  =  [' //ch_out  (1)  ( :nc_out (1) ) // 
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endif 


',V/ch_out(2)  (:nc_out  (2) ) // 

msecs;  Min  &  Max  Frequency  =  ' // 
ch_out (3) ( :nc_out (3) ) //'  &  '// 
ch_out (4)  ( :nc  out  (4))//'  Hz.' 


Check  for  finite  or  semi-infinite  line 

if (isemi.eq.O)  then 
Finite  length 

if (screen)  write (*, '  (/a/) ' ) 

+  USING  FIRST  GAUGE  PROPERTIES  TO  DEFINE  NORMALIZED  ANALOGUE  ' 
rhocd  =  rck(l)  /  xkod(l) 

Set  gamma  determining  forcing  function 
ffunc  =  twopi  *  fmax  *  rhocd  /  xkod(l) 
ffunc  =  dsqrt (ffunc) 
else 

Semi-infinite  length 
if (screen)  then 
numerr  =  0 

3  fmin  =  xkod(l)**2  /  (rck(l)  *  twopi) 

write  (*,'  (a, lpell.4,a, $) ' ) 

+  'ENTER  MINIMUM  FREQUENCY  (HZ)  FOR  ANALOGUE  [DEF='  fmin  '1  ' 
read(*,' (a)' ,end=2O0)  fpos  >  \ 

if(fpos.ne.'  '  .and.  fpos.ne.'/')  then 
rewind  19 

write  (19,'  (a)')  fpos ( :lnblnk  (fpos) ) 
rewind  19 

read(19, *, err=120)  fmin 
endif 

Check  for  legal  frequency  limits 
if (fmin. ge. fmax)  then 

write (*,'  (a)')  'FMIN  MUST  BE  LESS  THAN  FMAX.  TRY  AGAIN  ' 
go  to  120 
endif 

if (fmin.le.OdO)  then 

write (*,'  (a)')  'FMIN  MUST  BE  GP^ATER  THAN  ZERO.  TRY  AGAIN 
go  to  120  •n.vaAXiM. 

endif 

endif 

Set  ganuna  determining  forcing  function 
ffunc  =  dsqrt (fmax /fmin) 
endif 


Now  for  both  get  required  number  of  stages 

125  if  (screen)  then 
numerr  =  0 
130  nst  =  9 

write  (*,  '  (a,i2,a,$) '  ) 

+  'ENTER  REQUIRED  NUMBER  OF  STAGES  [DEF='  nst  '1  ' 

read(*, ' (a) ' ,end=200)  fpos 
if (fpos.ne.'  '  .and.  fpos.ne.'/')  then 
rewind  19 

write  (19, '  (a)' )  fpos (: Inblnk (fpos) ) 
rewind  19 
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read(19, *, err=130)  nst 
endif 

if (nst . gt .maxst )  then 
write (*,' (/2(a,i5),a)') 

+  'NUMBER  OF  STAGES  IS  ',nst,'.  MAXIMUM  IS  maxst  '  ' 
go  to  130 
endif 
else 

nst  =  9 
endif 


Solve  iteratively  for  gamma 


gam  = 
dgaiti  = 
do  while 
gamn 
gamn 
dgam 
gam 
enddo 


1.4d0 

IdO 

(dgam. gt . ld-10) 

—  IdO  —  ffunc  *  (IdO  —  gam)  /  (IdO  +  dsqrt (gam)  ) 
=  gamn** (IdO/nst) 

=  dabs (gam- gamn) 

=  gamn 


^®ta  —  IdO  /  (IdO  +  dsqrt (gam) ) 


if (screen)  then 
numerr  =  0 

Inform  user  of  gamma  value  and  ask  if  okay 
140  write(*, '  (a,  lpel3.6,a) ' )  ' CAPACITOR  GEOMETRIC  VALUE  IS  '  ,  gam, 

acceptable  (Y/N)?  [DEF=Y]',. 

^^(i*sq.78  .or.  i.eq.llO)  go  to  125  !  N  or  n 

if(i.ne.l3  .and.  i.ne.89  .and.  i.ne.l21)  go  to  140 
endif 


The  number  of  stages  and  their  distribution  is  now  known. 

Set  network  non-dimensionalize  by  the  first  resistor  and 
capacitor  values.  These  real  values  will  be  introduced  when 
the  governing  equations  are  solved. 

Note:  This  structure  assumes  that  the  same  "spacing"  of 

elements  will  be  used  for  all  cases  run,  and  thickness  of 
the  kapton/property  variations  will  be  introduced  when 
each  of  the  cases  to  be  run  is  considered  (i.e.  through 
R1  and  Cl  since  total  resistance  =  D/K,  and 
total  capacitance  =  RHO*C*D) 


......  ,  ,  .  - -*-*-.«-*.* A  O  JL/  J  • 

**************************i,-)r*******-k*i,*ir******-k*-h****-)r*-Hc-lri,***^r** 


* 

★ 

★ 


Now  set  the  first  capacitor/resistor  values  to  1.  The  network 
equations  will  become  non-dimensional  to  R1*C1,  and  R1*C1  will  be* 
^^troduced  when  the  time  steps  are  done . 

cap(l)  =  IDO 
res(l)  =  IDO 


^ssistance/capacitance  values  with  R/Rl  and  C/Cl 
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* 

'k 

* 

* 

* 

★ 


* 


Note,  in  real  terms: 

C(I)  =  GAM*C(I-1) 

R(l)  =  (R/C)  *C1/ (1+SQRT  (GAM)  ) 

R(I)  =  (R/C) *CI/SQRT (GAM) 

R(NST+1)  =  (R/C) *CNST*SQRT(GAM)/( 1+SQRT (GAM)) 
R(NST+1)  =  infinite  for  semi-infinite 


finite  length 


do  i=2,nst 

cap(i)  =  gam  *  cap(i-l) 

res(i)  =  dsqrt (gam)  *  (IdO  +  dsqrt(gam))  *  cap(i-l) 
enddo 

res(nst+l)  =  dsqrt (gam)  *cap (nst) 
if (isemi.eq.l)  res(nst+l)  =  ld38 


Get  total  resistance  and  capacitance  (normalized  by  R1,C1) 

rtot  =  ODO 
ctot  =  ODO 
do  i=l,nst 

rtot  =  rtot  +  res(i) 
ctot  =  ctot  +  cap(i) 
enddo 

if (isemi.eq.O)  rtot  =  rtot  +  res(nst+l) 


****************-k****ir*-k******-^-ir*-t;-k*-t!***-k-lt*******ifk*ic 
Now  fill  the  B-matrix  (i.e.  D  and  E  matrices) 

******  *********i,*******-k**i,ic********i,***ir**i,ic******i,i. 


Note:  Bi, i  =  -  1/Ci  *  (1/Ri  +  1/Ri+l) 

Bi, i+1  =  Bi+l,i  =  1/ (Ri+1*SQRT (Ci*Ci+l) ) 

do  i=l,nst 

d(i)  =  -(IdO  /  res (i)  +  IdO  /  res (i+1))  /  cap(i) 

if (i. It. nst)  e(i+l)  =  IdO  /  (res (i+1) *dsqrt (cap (i)  *cap (i+1) ) ) 


Set  input  transformation  matrix  A  to  identity  matrix 
(On  return  from  VALVAR  it  contains  the  eigenvectors.) 

do  j=l,maxst 
do  i=l,maxst 
a(i,j)  =  OdO 
enddo 

a(j,  j)  =  IdO 
enddo 

Now  get  eigenvector  matrix.  Eigenvalues  are  stored  in  matrix  D 
hence,  D  is  destroyed. 

call  valvar (d, e, a,  ic,  nst, maxst) 

if(ic.ne.O)  then 
numerr  =  9 
if (screen) 

go^to^20C)  ^  'Convergence  error  in  VALVAR  ~  stopping.' 


75 


Nov  5  09:55  1993  acq.F  Page  9 


endif 


************* ***********************^^^*^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 
This  completes  the  section  on  setting  up  the  non-dimensional 
analogue  network*  Now  we  get  into  dimensional  numbers  and 
*  solve  the  equations. 

**************************iri,*********->r****ir*******ir************i,**^^. 


***************** **********************^*^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 
Get  physical  properties  so  that  analogue  really  simulates  the 
*  heat  transfer  environment. 

**********************************i,***i,****************************** 


Get  RHO  C*K  If  semi-infinite  or  K/-D  and  RHO*C*D  if  finite  length, 
(i.e.  Input  voltage  =  temperature  =>  input  current  =  heat 
transfer  rate.) 


if (screen)  then 
numerr  =  0 
if (iserai.eq.O)  then 
First  get  start  time 

150  write(*,' (a,$)' )  'ENTER  START  TIME  (MS)  =>  ' 
read(*,*,err=150,end=200)  st 

St  =  St  /  ld3 

else 

St  =  OdO 

write(*, ' (a) ' )  'START  TIME  =  0  MSECS.' 
endif 

Get  end  time 

160  write(*, ' (a, $) ' )  'ENTER  END  TIME  (MS)  =>  ' 

read(*,*,err=160,end=200)  ft 

ft  =  ft  /  ld3 

Get  screen  print  interval 

170  write (*, ' (a, $) ' )  'SCREEN  OUTPUT  EVERY  N  STEPS,  ENTER  N  =>  ' 
read(*,*,err=170,end=200)  ipr 
if(ipr.le.O)  ipr  =  1 
write (*,  '(/)') 
endif 


Read  through  the  files  and  get  data 
do  nnn=l,npos 

Set  driving  voltages  to  zero 

vO  =  OdO 

vl  =  OdO 

vOb  =  OdO 

vlb  =  OdO 

v5b  =  OdO 

vforeb  =  OdO 

vforlb  =  OdO 

vfor5b  =  OdO 

if (npos .gt . 1)  then 

Set  temperature  data  file  names 
fname  =  name(nnn) 
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ncstem  =  Inblnk  (fnaitie) 
foutt  =  fname ( rncstem) //' T. 1' 
foutb  =  fname (: ncstem) //'B.l' 
foutq  =  fname (: ncstem) //'Q. 2' 
endif 

if (screen)  then 

Tell  user  the  files  we  are  reading  and  writing 
write(*,' (a)')  'READING  FILE:-  ' //foutt (1 : length (foutt ) ) 
if (isemi.eq.O)  write (*, '  (a) '  )  ' READING  FILE : -  ' // 

foutb ( 1 : length ( foutb ) ) 

write(*,'  (a)')  '  WRITING  FILE:-  ' //foutq(l:length(foutq) ) 
endif 

Open  upper  surface  temperature  file  for  read 
open (1, file=foutt, status=' old' , form=' unformatted' , 
recl=1024 , access='  direct' ) 

read (1,  rec=l)  trdr 
read(l,  rec=2)  tcdr 

Compute  required  properties 
if  (isemi.eq.O)  then 
Finite-length 
rhocd=rck (nnn) /xkod (nnn) 
rl  =  IdO/ (xkod (nnn) *rtot) 
cl  =  rhocd/ctot 
else 

Semi-infinite  length 

Compute  first  stage  capacitance  and  resistance 
cl  =  IdO/dsqrt (twopi*fmax*beta*beta/rck (nnn) ) 
rl  =  IdO/ (twopi*fmax*beta*cl) 
endif 

rrlcl=ldO/ (rl*cl) 
ichnum  =  ichar <tcdr.iss) 

if (ichnum.gt . 96)  tcdr.lss  =  char (ichnum-32) 
ntimes  =  tcdr.nsampl 
if  (ipvalt.gt.O)  call  free(ipvalt) 
ipvalt  =  malloc (ntimes*4) 
ntrecs  =  2  +  (ntimes  +  255)  /  256 

j2  =0 

do  irec=3, ntrecs 
jl  =  j2  +  1 

j2  =  minO ( j2+25 6,  ntimes ) 

read(l,  rec=irec)  (vt ( j) , j=jl, j2) 
enddo 

close (1,  status=' keep' ) 
if (iptime . gt . 0 )  call  free (iptime) 
iptime  =  malloc (ntimes*4) 
if (tcdr.lss. eq. 'I' )  then 

inquire (file=' INDEXER. 1' , exist=exists) 
if  (exists)  then 

open (1, file='  INDEXER. 1'  , status='  old' ,  access=' direct' , 
form=' unformatted' ,recl=1024) 
nirecs  =  2  +  (ntimes  +  127)  /  128 
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j2  =0 

jO  =  0 

jn  =0 

do  irec=3, nirecs 
jl  =  j2  +  1 

j2  =  minO ( j2+128, ntimes ) 

ninput  =  j2  -  jl  +  1 

read(l,rec=irec)  (r8time(j) , j=l, ninput) 
k  =0 

do  j=jl,j2 

k  =  k  +  1 

rt(j)  =  r8time(k) 
if  (rt ( j) .le.st)  jO  =  j 

if  (rt  (j)  .le.ft)  jn  ?=  j 

enddo 
enddo 

close (1, status=' keep' ) 
else 

numerr  =  3 

if(screen)  write (* , ' (a) ' ) 

'FILE  " INDEXER. 1"  DOES  NOT  EXIST  —  STOPPING.' 

go  to  200 
endif 

elseif (tcdr.lss.eq.'C' )  then 
deltat  =  IdO  /  tcdr.lsr 
jO  =1.5+  St/  deltat 

jii  =  1.5  +  ft  /  deltat 
jO  =  maxO (1, jO) 

jn  =  minO (ntimes, jn) 
do  j=l,  ntimes 

rt(j)  =  (j  -  1)  *  deltat 
enddo 
endif 


Check  for  meaningful  time  interval  indices 
numerr  =  4 

if(j0.1t.l  .or.  jO.gt. ntimes  .or.  jO.ge.jn)  go  to  200 
if(jn.lt.l  .or.  jn.gt .ntimes)  go  to  200 
nud  =  jn  -  jO  +  1 

nub  =  jn 

Open  heat  transfer  rate  file  for  write 

open (3,  file— foutq, status=' unknown' , form=' unformatted' , 
recl=1024 ,  access=' direct ' ) 
qrdr  =  trdr 

qcdr  =  tcdr 

qcdr.unityp  =  5 
qcdr. unilab  =  'Watts/m''2' 
qcdr.ctype  =  2 
nclabl  =  Inblnk (tcdr. label) 
qcdr. label (nclabl:nclabl)  =  ' Q' 
nstory  =  minO (Inblnk (qcdr . story) , 66) 
qcdr. story  =  qcdr. story ( rnstory) //'  QDOT  FROM  ACQ' 
write (3, rec=l)  qrdr 
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write ( 3, rec=2)  qcdr 
if (ipvalo.gt.O)  call  free(ipvalo) 
ipvalo  =  malloc (ntimes*4) 

if (isemi.eq. 0)  then 

Open  lower  temperature  file  for  read 
open  (2,  file=foutb,  status='  old' ,  form='  unformatted' , 
recl=1024, access=' direct' ) 
read (2, rec=l)  brdr 
read (2, rec=2)  bcdr 
ichnum  =  ichar (bcdr. Iss) 

if (ichnum.gt.96)  bcdr. Iss  =  char (ichnum-32) 

if (bcdr . Iss . ne . tcdr . Iss  .or.  bcdr.nsampl.ne.tcdr.nsampl)  then 
numerr  =  5 

if (screen)  write (*, ' (a) ' ) 

'TOP /BOTTOM  TEMPERATURE  FILE  INCOMPATIBILITY  —  STOPPING  ' 
go  to  200 
endif 

if (ipvalb.gt .0)  call  free (ipvalb) 
ipvalb  =  malloc (ntimes*4) 
j2  =0 

do  irec=3 , ntrecs 
jl  =  j2  +  1 

j2  =  minO ( j2+2S6, ntimes) 

read (2, rec=irec)  (vtb( j) , j=jl, j2) 
enddo 

close (2, status=' keep' ) 
endif 


★ 


*■ 


★ 


Set  boundary  conditions  based  on  length  analogue  chosen 


if (isemi.eq. 0)  then 
Finite 
sumres  =  OdO 
do  m=l,nst 


sumres  =  sximres  +  res  (m) 

t (m)  —  vt(jO)  +  (vtb(jO)  —  vt(jO))  *  rl  *  xkod(nnn) 
enddo 

do  m=l,nst 
y(m)  =  OdO 
do  i=l,nst 

y (m)  =  y (m)  +  a(i,m)  *  dsqrt (cap (i) )  *  t(i) 
enddo 
enddo 
else 


Semi-inf inite 
do  m=l,nst 
y(m)  =  OdO 
do  i=l,nst 


y  (m) 

=  y(m) 

enddo 

y  (m)  = 

vt  (jO) 

enddo 

vOb  = 

vt ( j  0 ) 

vlb 

vOb 

v5b 

.5d0  * 

+  a{i,m)  *  dsqrt (cap (i) ) 
*  y(m) 


(vOb  +  vlb) 


*  sumres 
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vforeb  =  vOb/ (res (nst+1) *dsqrt (cap (nst ) ) ) 
vforlb  =  vlb/ (res (nst+1) *dsqrt (cap (nst) ) ) 
vfor5b  =  v5b/ (res (nst+1) *dsqrt (cap (nst) ) ) 
endif 

Set  up  loop  counter  (lOUT  =  output  file  counter) 
iout  =  0 

Loop  for  picking  up  time  and  values 

(use  upper  surface  temperature  file  as  timing  file) 

if(nud.eq.O)  then 

Inform  user  there  are  no  data 
if (screen)  write (* , '  (a) '  ) 

+  'NO  TEMPERATURE  VALUES  IN  REQUESTED  TIME  INTERVAL.' 
else 

ul  =  OdO 

do  i=l,nst 

ul  =  ul  +  a(l,i)  *  y(i) 

enddo 

vto(jO)  =  (vt(jO)  -  ul  /  dsqrt (cap (1) ) )  /  rl 
do  iq=j0+l, jn 
time  =  rt(iq) 

vO  =  vt (iq-1) 

vl  =  vt (iq) 

Get  time  step 

dt  =  rt(iq)  -  rt(iq-l) 

If  finite  length  line,  get  lower  temp  at  this  and  last  time 
if (isemi.eq. 0)  then 
vOb  =  vtb(iq-l) 
vlb  =  vtb (iq) 
endif 

Now  we  have  time  and  value (s)  so  solve  for  Qdot 
Assume  that  voltage  varies  linearly  between  samples 

*****-k*****ir*-k***i!*********-k**it*ic******ir*-k-kic-k*ir**ie* 

v5  =  .5d0  *  (vO  +  vl) 

vforce  =  vO/ (res (1) *dsqrt (cap(l) ) ) 
vforcl  =  vl/ (res (1) *dsqrt (cap(l) ) ) 
vforc5  =  v5/ (res (1) *dsqrt (cap (1) ) ) 
if (isemi.eq. 0)  then 

v5b  =  ,5d0  *  (vOb  +  vlb) 
vforeb  =  vOb/ (res (nst+1) *dsqrt (cap (nst) ) ) 
vforlb  =  vlb/ (res (nst+1) *dsqrt (cap (nst) ) ) 
vfor5b  =  v5b/ (res (nst+1) *dsqrt (cap (nst) ) ) 
endif 

Now  march  through  all  the  nodes  to  get  solution  at  time  n+1 
do  i=l,nst 

Use  a  fourth-order  Runge-Kutta  operation 
dy(I)/dt  =  Eigenvalue (I) *y (I)  +  Eigenvector (1, I) *VFORCE 
+  Eigenvector (N, I) *VFOREB 

=  f(t,y)  (If  finite  length  line) 

First  step:  An  =  dt*f(t,y) 
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★ 

★ 

★ 

★ 


Second  step:  Bn  =  dt*f {t+dt/2,  y+An/2) 

Third  step:  Cn  =  dt*f (t+dt/2, y+Bn/2) 

Fourth  step:  Dn  =  dt*f (t+dt , y+Cn) 

Last  step:  y(t+dt)  =  y(t)  +  (An  +  2Bn  +  2Cn 

y(i)  +  a(l,i)  * 


an 

=  rrlcl  *  dt  *  (d(i) 

+ 

+  a  (nst, i)  *  vforeb) 

bn 

=  rrlcl  *  dt  *  (d(i) 

+ 

+  a(l,i)  *  vforc5  + 

cn 

=  rrlcl  *  dt  *  (d(i) 

+ 

+  a(l,i)  *  vforcS  + 

dn 

=  rrlcl  *  dt  *  (d(i) 

+ 

+  a(nst,i)  *  vforlb) 

y(i) 

=  y(i)  +  (an  +  2d0  *  ( 

+  Dn)/6 
vforce 


*  (y(i) 
a (nst, i) 

*  (y(i) 
a (nst, i) 


+  .5d0  *  an) 

*  vforSb) 
+.5d0  *  bn) 

*  vforSb) 

*  (y(i)+cn)  +  a(l,i) 

m  +  cn)  +  dn)  /  6d0 


vforcl 


enddo 


Now  transform  back  to  get  voltage  at  station  1 

ul  =  OdO 

do  i=l,nst 

ul  =  ul  +  a(l,i)  *  y(i) 

enddo 

voll  =  ul  /  dsqrt (cap (1) ) 

And  to  current  (i.e.  heat  transfer  rate) 
amp  =  (vl  -  voll)  /  rl 

Prepare  output  file  array 
vto(iq)  =  amp 
if (screen)  then 
iout  =  iout  +  1 

if (mod (iout, ipr) . eq, 0)  write (*, '  (2 (a, lpel4 . 6) , a) '  ) 
+  'Time  =',time,'  secs;  Heat-f lux=' , amp, '  W/sq  m.' 

endif 

Loop  back  for  next  value 
enddo 


do  j=l, jO-1 
vto  ( j )  =  0.0 
enddo 

do  j= jn+1, ntimes 
vto(j)  =  0.0 
enddo 

j2  =  0 

do  irec=3, ntrecs 
jl  =  j2  +  1 

j2  =  minO ( j2+256, ntimes) 

write (3, rec=irec)  (vto(j),j=jl,j2) 
enddo 

close  (3, status=' keep'  ) 

if ( .not . screen)  then 
Write  to  Log  file 
if (isemi . eq. 0)  then 

write (4, '  (a) ' )  '*  Heat  flux  ' //foutq ( : Inblnk (foutq)  ) // 
'  from  ' //foutt (: Inblnk (foutt) ) // 

'  and  '/ /foutb (: Inblnk  (foutb)) 

else 
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write  (4, '  (a) ' )  '*  Heat  flux  ' //foutq( :lnblnk (foutq) ) // 
'  from  ' //foutt (:lnblnk(foutt) ) 

endif 

endif 

endif 

Loop  back  for  next  file 
enddo 
Fini 

i  =  ieee_flags ("clear", "exception", "all",  fpos) 

call  exit (numerr) 
end 
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